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ABSTRACT 

We revisit the classical problem of synthesizing spectral properties of a galaxy using 
a base of star clusters, approaching it from a probabilistic perspective. The problem 
consists of estimating the population vector x, composed by the contributions of n^, 
different base elements to the integrated spectrum of a galaxy, and the extinction Ay , 
given a set of absorption line equivalent widths and continuum colors. The formalism 
is applied to the base of 12 elements defined by Schmidt etal. (1991) as correspond- 
ing to the principal components of the original base employed by Bica (1988), and 
subsequently used in several studies of the stellar populations of galaxies. The explo- 
ration of the 13-D parameter space is carried out with a Markov chain Monte Carlo 
sampling scheme based on the Metropolis algorithm. This produces a smoother and 
more efficient mapping of the P(x, Ay) probability distribution than the traditionally 
employed uniform-grid sampling. 

This new version of Empirical Population Synthesis is used to investigate the abil- 
ity to recover the detailed history of star formation and chemical evolution using this 
spectral base. This is studied as a function of (1) the magnitude of the measurement 
errors and (2) the set of observables used in the synthesis. Extensive simulations with 
test galaxies are used for this purpose. Emphasis is put on the comparison of input 
parameters and the mean x and Ay associated with the P{x,Av) distribution. It is 
found that only for extremely low errors {S/N > 300 at 5870 A) aU 12 base propor- 
tions can be accurately recovered, though the observables are recovered very precisely 
for any S/N. Furthermore, the individual x components are biased in the sense that 
components which carry a large fraction of the light tend to share their contribution 
preferably among components of same age. Old, metal poor components can also be 
confused with younger, metal rich components due to the age x metallicity degener- 
acy. These compensation effects are linked to noise-induced linear dependences in the 
base, which very effectively redistribute the likelihood in x-space. The age distribution, 
however, can be satisfactorily recovered for realistic S/N (~ 30). We also find that 
synthesizing equivalent widths and colors produces better focused results that those 
obtained synthesizing only equivalent widths, despite the inclusion of the extinction 
as an extra parameter. 
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1 INTRODUCTION 

The stellar populations of galaxies carry a record of their 
star-forming and chemical histories, from the epoch of for- 
mation to the present. Their study thus provide a powerful 
tool to explore the physics of galaxy formation and evolu- 
tion. Despite the detailed understanding of stellar evolution 
and of the properties of single stellar populations such as 
star clusters achieved over the past century, uncovering the 
history of galaxies through the study of their integrated light 
has proven to bo a difficult problem (Pickles 1985; Worthey 
1994 and references therein). 

Two global approaches have been developed to tackle 
this issue. The first of these is Evolutionary Population Syn- 
thesis, which performs ab initio calculations of the spectral 
evolution of gala^xies on the basis of stellar evolution the- 
ory, stellar spectral libraries, and prescriptions for the ini- 
tial mass function, star formation rate and chemical evolu- 
tion (Tinsley 1972; Larson & Tinsley 1974, 1978; Guider- 
doni & Rocca-Volmcrange 1987; Fioc & Rocca-Volmerange 
1997; Chariot & Bruzual 1991; Bruzual & Chariot 1993; Lei- 
therer etal. 1996, 1999 and references therein). The other 
approach is Empirical Population Synthesis (EPS here- 
after), also known as 'Stellar Population Synthesis with a 
Database', which uses observed properties of stars or star 
clusters as a base on which to decompose the population mix 
in galaxy spectra (Spinrad & Taylor 1971; O'Connell 1976; 
Faber 1972; Pickles 1985; Bica 1988; Pelat 1997, 1998; Bois- 
son et al. 2000) . Of particular interest is the EPS method of 
synthesizing the spectra of galaxies using a base of star clus- 
ters. Perhaps the most appealing property of this method is 
the fact that the results do not hinge on assumptions about 
stellar evolution tracks nor initial mass functions, as these 
are 'mother nature given' in the empirically built base spec- 
tra. This technique has been pioneered by the work of Bica 
(1988, hereafter B88), where the stellar content of a large 
sample of galaxies was synthesized using the spectral base 
of star clusters described in Bica & AUoin (1986a, 1986b, 
1987). Since then, several studies made use of this technique 
(Schmidt, Bica & Dottori 1989; Bica, Alloin & Schmidt 1990; 
Jablonka, Alloin & Bica 1990; Schmidt, Bica & Alloin 1990; 
Jablonka & Alloin 1995; De Mello et al. 1995; Bonatto et al. 
1996, 1998, 1999, 2000; Schmitt, Storchi-Bergmann & Cid 
Fernandes 1999; Kong & Cheng 1999; Raimann etal. 2000). 

In its original version (B88), this method made use of 
a large base of n* = 35 elements, combined in different pro- 
portions (.Tj: 1 = 1... 35) to synthesize a number of equiv- 
alent widths (W) of conspicuous absorption features. An 
uniform sweep of the parameter space along paths consis- 
tent with simple predefined chemical evolution scenarios was 
performed, and all combinations which produced W's within 
10% of the observed values were considered 'solutions'. The 
final population vector x was talcen to be the arithmetic 
mean of all such sampled solutions. 

One of the main concerns raised by B88 method is 
that of the uniqueness of the solution, as first discussed by 
Schmidt etal. (1991, hereafter S91). With a 35-D popula- 
tion vector and at most 9 W's to be synthesized, one has a 
highly degenerate algebraic system. As shown by S91, how- 
ever, the original base was highly redundant, with several 
linearly dependent elements, and others which were practi- 
cally so, as deduced by a principal component analysis. S91 
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then proposed the use of a reduced base of 12 elements. S91 
also criticized B88's uniform sampling of the x space (a 'dis- 
crete combination procedure' in their terms) and the use of a 
'mean solution'. Instead, they developed a constrained min- 
imization procedure, which, as shown by their simulations 
with test galaxies built out of the base, is able to retrieve 
with good accuracy the correct population vector. Another 
difference of the S91 method compared to B88 is the possi- 
bility to search for solutions along the entire age x metal- 
licity plane. B88 assumed that the (iiciiiical enrichment of 
the galaxy occurred during its formation (-^ 10 Gyr ago), 
with all stars younger than that having the same Z as the 
oldest most metal rich component. Although this may be a 
reasonable approximation for the nuclear region of isolated 
galaxies, it does not allow for the effects of mergers, inflows 
or outflows, and is not applicable to high redshift sources, 
where the spectra usually integrate over the entire galaxy, 
including regions with different star formation and chemical 
histories (Jablonka etal. 1990; Jablonka & Alloin 1995). 

The next improvement in this technique was the intro- 
duction of continuum colors in the synthesis. This was done 
by Schmitt, Bica & Pastoriza (1996), who, as B88, use a 
discrete combination procedure search for acceptable solu- 
tions, sweeping the parameter space along chemical enrich- 
ment paths. This version of EPS has been used in several 
subsequent studies (Bonatto etal. 1996, 1998, 1999, 2000; 
Kong & Cheng 1999). Another step was taken by Bonatto, 
Bica & Alloin (1995), who extended the star cluster base to 
the space ultraviolet (1200-3200 A). They too use an arith- 
metic average to represent the solution of the synthesis. 

In two recent papers, Pelat revisited the EPS problem 
of synthesizing H^'s using an elegant formalism based on 
convex algebra concepts. In Pelat (1998) he shows how the 
algebraic solution of the EPS problem can be narrowed down 
to a sub-space of x delimited by a set of 'extreme-solutions', 
any convex combination of which yields an exact solution 
in the underdetermined cases (more parameters than con- 
straints). His results shed new light into the issue of alge- 
braic degeneracy, which has always haunted EPS due to the 
widespread belief that with more parameters than observ- 
ables one is bound to fit the data in one way or another. 
He proposes a quick test of whether the VF's of a galaxy are 
synthesizable by computing the region in VF-space spanned 
by all physical combinations of the base elements. We veri- 
fied that in the case of the 12 elements base used by S91 this 
region is narrow (Leao 2001), such that small measurement 
errors are enough to place objects outside this 'synthetic 
domain', hence preventing the existence of mathematically 
exact solutions. In real applications, one therefore expects 
the degeneracy of EPS to be more of a statistical nature 
than algebraic. In these cases, as well as in overdetermined 
problems, Pelat (1997) demonstrates that the best model is 
to be found along the boundary of the synthetic domain by 
minimizing an appropriately defined x^-liko distance. Un- 
certainties in the population vector can also be readily com- 
puted, as shown by Moultaka & Pelat (2000). 

Boisson et al. (2000) recently applied Pelat's method to 
a sample of 12 active galaxies from Serote Roos et al. (1998). 
They synthesize nw = 47 absorption lines with a base of 
n« = 30 stars from the stellar spectral libraries of Serote 
Roos, Boisson & Joly (1996) and Silva & Cornell (1992). 
Their results nicely illustrate the power of this new EPS 
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technique, which will certainly be extended to larger sam- 
ples and other classes of galaxies in future studies. The B88 
method and its variants, on the other hand, already have a 
large body of literature associated to it. Nevertheless, this 
popular and attractive EPS technique, while intuitively ac- 
ceptable, has neither been adequately tested nor cast onto a 
sound mathematical formalism which supports its applica- 
tion. Furthermore, it gives no measure of the uncertainties 
or potential biases on the population vector, thus giving no 
simple way to prevent overinterpretations of the resulting 
population of the age x metallicity plane. These issues are 
the focus of the present paper. 
Our main goals are to: 

(i) Elaborate a mathematically consistent formulation of 
the EPS problem in the context of probability theory. 

(ii) Develop and test a sampling method to explore the pa- 
rameter space in EPS problems. In particular, we aim to 
improve upon the Schmitt etal. (1996) technique of syn- 
thesizing both equivalent widths and colors, but within the 
context of a probabilistic formulation. 

(iii) Use the method to investigate the effects of measure- 
ment errors and of the specific set of observables used upon 
the results of the synthesis with the reduced B88 star clus- 
ter base described in S91. In particular, we wish to evaluate 
whether these factors introduce biases in the population mix 
inferred from this popular synthesis technique and at what 
level of detail can its results be trusted. 

This paper is organized as follows: Item fi) of the above 
list is addressed in section |^. In section ^ we deal with 
point (ii). This is done by investigating the applicability of a 
Metropolis algorithm to sample the parameters probability 
distribution in EPS, and testing the method with simula- 
tions based on Bica's spectral base. In section ^ we address 
item (iii) by means of extensive simulations which explore 
the ability to recover the detailed star-formation and chem- 
ical histories of galaxies with this base. Finally, section ^ 
summarizes our main results. 



2 FORMALISM 
2.1 Basic equations 

EPS studies seek to find combinations of a spectral base 
which reproduce a given set of measured observables, of- 
ten taken as the equivalent widths Wj of nw conspicuous 
absorption features. The base consists of elements repre- 
senting well defined simple stellar populations such as star- 
clusters, with equivalent widths W*j and corresponding con- 
tinuum fluxes over the lines F*^ (j = 1, . . . nw, i = 1, . . . n*) 
normalized at a reference wavelength. Denoting by Xi the 
fractional contribution of the i-th base element to the total 
flux at the reference wavelength, one obtains a system of nw 
equations 



, nw 



(1) 



for the synthetic W's (e.g., Joly 1974). This set of constraints 
can be augmented by synthesizing nc observed continuum 
fluxes Ck {k = l,...nc), provided allowance is made for 
reddening, here parametrized by the V-band extinction Ay- 



Assuming all populations are equally reddened, one ob- 
tains 

Ck^Cki:>i,Av)=gk{Av)^C*kXi; fc = l,...,nc (2) 



where a distinction is made between C* and F* , since the 
nc fluxes to be synthesized need not correspond to the wave- 
lengths of the nw absorption lines. We hereafter refer to the 
Cfc's as continuum colors, as they are in fact ratios of contin- 
uum fluxes with respect to the continuum at the normaliza- 
tion wavelength. The function gk{Av) reddens the normal- 
ized color at wavelength by A\^ according to a specifled 
reddening law. 

Finally, physical solutions must further satisfy the nor- 
malization and positivity constraints: 



E 



Xi 



1; Xi > for i — 1, n*, and Av > 0. 



(3) 



The normalization condition effectively reduces one de- 
gree of freedom. When modeling colors, however, one has to 
introduce Av as a further parameter, so that the number of 
parameters still is n*. 

2.2 Probabilistic formulation 

The data V to be modeled are thus composed of a set of 
nobs ~ nw + nc observables. The measurement errors in 
each observable, collectively denoted by a, are assumed to 
be known. Given these, the problem of EPS is to estimate 
the population vector x and the extinction Av that 'best' 
represents the data according to a well deflned probabilistic 
model. It is equally important to estimate the uncertainties 
in the parameters, as these prevent over-interpreting the re- 
sulting mixture of stellar populations at a level of detail not 
warranted by the data or by intrinsic limitations of the base. 

The probability of a solution (x, Av) given the data D 
and the errors a, is given by Bayes theorem: 



Pix,Av\V, a,H) 



P(©|x, Av,H)P{y., Av\a, H) 



(4) 



p{v\n) 

Ti summarizes the set of assumptions on which the in- 
ference is to be made. They include: the mapping between 
parameters and observables is given by eqs. (hi) and (eI); x 
and Av must satisfy the constraints expressed in eq. rtsp; the 
stellar population in the target galaxy is well represented by 
the base elements; the observational errors are Gaussian; the 
reddening law is known. 

The likelihood P(I'|x, Av, o", H) is a measure of how 
good (or bad) is the fitting to the data for model parameters 
X and Av- Under the hypothesis of Gaussian errors, the 
probability of the data given the parameters is 



P{V\x,Av,(T,n) oc e~^, 

with £ defined as half the value of 



£{x,Av) = 



ix^(x,A.)4^ 



H/,-(x) 



+ 



nc 



3 = 1 

Cf--a(x,Av) 



(5) 



(6) 
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The likelihood, as the total x^, separates into a W and a 
color related term, of which only the latter depends on Ay- 

The normalizing constant P(X'|7i) in Bayes theorem 
(the 'evidence of H') is irrelevant to the level of inference 
discussed here, i.e., the estimation of x and Ay- It can be 
used, for instance, to compare different spectral bases. 

P(x, v4\^|7^) is the joint a priori probability distribution 
of X and Ay, and states what values the model parame- 
ters might plausibly take. For instance, physically accept- 
able population fractions should have priors that are zero 
in regions of the parameter space where the constraints of 
positivity and unit sum are not satisfied. A prior that is uni- 
form on the parameters and includes these constraints is a 
non-informative prior, because it does not impose any con- 
straints on the solutions besides those expressed by eq. 
This is the prior that will be used in this work. 

For a non-informative prior the posterior probability 
P(x, j4\/|©, a, Ti.) is simply proportional to the likelihood: 



P{x,Av\V,a,H) oc e 



-£(x,Av) 



(7) 



This expression contains the full solution of the EPS 
problem, as embedded in it is not only the most proba- 
ble model parameters but also their full probability distri- 
butions. Furthermore, projected posterior distributions for 
any of the n* parameters can be obtained by marginalizing 
eq. ^ with respect to all other parameters. For instance, the 
probability density of proportion Xi of the i-th base element 



P{x,Av\'D,a,n)dxi 



dxi-idxi+i . . . dxn^dAv 



(8) 



and equivalently for Av- Similarly, one can construct joint 
posteriors for, say, the total proportion of all base compo- 
nents of same age or Z, or for the mean age and Z of the 
stellar population. This would of course provide a coarser 
description of a galaxy's stellar content than the individual 
posteriors P{xi\'D, a, H), but that may be all that is possible 
under some circumstances, such as when only a reduced set 
of observables is available or when observational errors are 
large (section ^ . 

It is easy to see how the B88 and Schmitt etal. (1996) 
synthesis techniques fit into this general probabilistic for- 
mulation. By synthesizing the observables within 10% 
'error boxes', these authors implicitly assumed a box-car 
likelihood function P(D|x, yly , cr, 7i), whereas by perform- 
ing arithmetic means over all accepted (x, Ay) combina- 
tions in a uniform grid they are implicitly sampling the cor- 
responding posterior probability distributions. Also, the a 
priori constraints on the occupation of the age x Z plane 
imposed by these authors (but relaxed in subsequent works) 
simply refiect their use of an informative prior. This gives 
some formal justification for their heuristically designed EPS 
method. It is therefore reasonable to expect that an imple- 
mentation of EPS based on the formalism presented here 
should give results roughly compatible with those previously 
obtained. An uniform sweep of the (x. Ay) space is however 
a very inefficient sampler for (^, so some improvement is 
needed there. This is discussed next. 



3 SAMPLING METHOD AND TESTS 

Despite its formal merits, at the computational level our 
probabilistic approach faces the same basic difficulty as pre- 
vious EPS codes, namely, the high dimensionality of the 
problem. For spectral bases with astrophysically interesting 
resolution in age and metallicity the number of elements 
n« quickly becomes large enough to render the exploration 
of the parameter space a non-trivial task. Before discussing 
methods to sample the (x. A y) s pace, we present the spec- 
tral base used in this work (j p.l| ) and desc ribe how we deal 
with the uncertainties in the observables (§3.2). 



3.1 The spectral base and observables 

The base used for the synthesis here is that of S91. It con- 
tains n* = 12 population groups, spanning five age bins — 10 
Gyr (which actually represent globular cluster-like popula- 
tions), 1 Gyr, 100 Myr, 10 Myr and HII (corresponding to 
current star formation, and represented by a pure Fx oc A^^ 
continuum based on the spectrum of 30 Dor) — and four 
metallicities— 0.01, 0.1, 1 and 4 Zq (Table 1). The observ- 
ables in this base comprise the nw — 9 equivalent widths of 
the absorption lines Call K A3933, CN A4200, G band A4301, 
Mgl A5175, Call A8543, Call A8662, US, H7 and H/3, as well 
as nc = 7 continuum fiuxes at selected pivot wavelengths: 
3290, 3660, 4020, 4510, 6630, 7520 and 8700 A, aU normal- 
ized to 5870 A. According to Bica & AUoin (1986a, 1987) 
and Bica, AUoin & Schmitt (1994), the equivalent width 
windows and continuum points were defined based on very 
high signal to noise spectra of galaxies, with the express goal 
of using them to synthesize the stellar population of galax- 
ies. These spectra were obtained by creating two average 
spectra, of blue and red galaxies, where the stellar popula- 
tion features can be clearly traced. The base values for these 
quantities have been previously published in Bica & AUoin 
(1986b, 1987) and Schmitt et al. (1996), and some were mea- 
sured from data in Bica etal. (1994). These are recompiled 
in Table 2 , as they bear a direct impact on the analysis 
below. There are thus a maximum of 9 -I- 7 = 16 observ- 
ables to model with = 12 parameters: n* — 1 population 
fractions plus Ay - However, in several of the tests presented 
below we shall make use of a reduced subset of observables, 
as observational data sets seldom cover the whole spectral 
range spanned by this base. 

According to Bica & AUoin (1986a, 1987) and Bica, Al- 
loin & Schmitt (1994), since the ultimate goal of the contin- 
uum points and W's is to synthesize the stellar population 
of galaxies, the W windows and continuum points were de- 
fined based on very high signal to noise spectra of galaxies. 
These spectra were obtained by creating two average spec- 
tra, of blue and red galaxies, where the stellar population 
featuress can be identified. 

The reddening of the colors is modeled with the extinc- 
tion law described in Cardelli, Clayton & Mathis (1989, with 
Ry = 3.1). 



3.2 Errors in the observables 

The measurement errors in the observables play a key role 
in determining the structure of the likelihood function. 
Whereas such errors are available when analyzing observed 
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Base Elements Used 



HII 10 Myr 


100 Myr 


1 Gyr 


10 Gyr 


log(Z/Z0) 


10 


8 


5 


1 


0.6 


12 11 


9 


6 


2 


0.0 






7 


3 


-1.0 








4 


-2.0 



Table 1. Ages, metallicities and numbering convention for the star clusters in the base. 



Equivalent Widths (A) 



# 


K 


ON 


G 


Mgl 


CaTi 


CaT2 


H<5 


H7 


H/3 


1 


21.1 


17.5 


11.1 


10.1 


6.8 


6.0 


4.4 


4.9 


3.5 


2 


17.3 


12.0 


9.3 


7.4 


5.5 


5.0 


4.4 


4.9 


3.5 


3 


11.0 


4.5 


5.9 


3.8 


3.4 


3.3 


4.4 


4.9 


3.5 


4 


4.7 


0.5 


2.5 


0.7 


1.2 


1.6 


4.4 


4.9 


3.5 


5 


17.3 


13.9 


9.0 


9.1 


6.8 


6.0 


9.7 


7.7 


7.5 


6 


14.0 


9.6 


7.4 


6.2 


5.5 


5.0 


9.7 


7.7 


7.5 


7 


8.9 


3.6 


4.6 


3.2 


3.4 


3.3 


9.7 


7.7 


7.5 


Q 
O 




O.VJ 


1.5 


3.2 


6.8 


6.0 


lU. 


Q Q 


7 Q 

/ .y 


9 


3.8 


2.2 


1.2 


2.5 


5.5 


5.0 


10.5 


9.9 


7.9 


10 


2.6 


1.4 


0.3 


2.5 


8.2 


6.9 


4.5 


3.5 


3.9 


11 


2.2 


1.1 


0.3 


2.0 


6.6 


5.8 


4.5 


3.5 


3.9 


12 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 








Continuum 


over the lines 








1 


U.o4 


U.4o 


0.55 


0.87 


1.03 


1.04 


U.44 


U.Oo 


U. (if 


2 


0.48 


0.59 


0.65 


0.90 


0.96 


0.96 


0.56 


0.67 


0.83 


Q 
O 


U. / z 


U. / 


0.81 


0.95 


0.84 


0.84 


U. ( 


U.oz 


U.yi 


4 


0.94 


0.96 


0.97 


1.01 


0.72 


0.71 


0.96 


0.97 


1.00 


tr 




i.Uo 


i.Uo 


1.06 


1.04 


0.70 


0.70 


i.Uo 


i.UD 


i.Uo 


6 


1.03 


1.05 


1.06 


1.04 


0.70 


0.70 


1.05 


1.06 


1.06 


7 


1.03 


1.05 


1.06 


1.04 


0.70 


0.70 


1.05 


1.06 


1.06 


8 


1.89 


1.75 


1.68 


1.23 


0.59 


0.58 


1.79 


1.65 


1.38 


9 


1.89 


1.75 


1.68 


1.23 


0.59 


0.58 


1.79 


1.65 


1.38 


10 


1.55 


1.34 


1.27 


1.06 


0.95 


0.94 


1.40 


1.25 


1.08 


11 


1.55 


1.34 


1.27 


1.06 


0.95 


0.94 


1.40 


1.25 


1.08 


12 


2.34 


2.03 


1.93 


1.37 


0.40 


0.38 


2.11 


1.90 


1.56 










Continuum Colors 










# 


3290 


3600 


4020 


4510 


6630 


7520 


8700 






1 


0.17 


0.27 


0.39 


0.71 


1.01 


1.01 


1.04 






2 


0.27 


0.37 


0.52 


0.77 


0.99 


0.96 


0.96 






3 


0.42 


0.52 


0.74 


0.88 


0.96 


0.89 


0.84 






4 


0.57 


0.67 


0.95 


0.99 


0.92 


0.81 


0.71 






5 


0.71 


0.65 


1.04 


1.08 


0.92 


0.81 


0.70 






6 


0.71 


0.65 


1.04 


1.08 


0.92 


0.81 


0.70 






7 


0.71 


0.65 


1.04 


1.08 


0.92 


0.81 


0.70 






8 


1.10 


0.94 


1.84 


1.52 


0.82 


0.68 


0.58 






9 


1.10 


0.94 


1.84 


1.52 


0.82 


0.68 


0.58 






10 


2.10 


1.52 


1.45 


1.11 


0.95 


0.99 


0.94 






11 


2.10 


1.52 


1.45 


1.11 


0.95 


0.99 


0.94 






12 


3.78 


2.56 


2.20 


1.73 


0.74 


0.56 


0.37 





Table 2. Observables in the base, which make the matrices W*-, F*. and C* . All continuum fluxes are normalized to the continuum at 
5870 A. 

galaxy spectra, they have to be postulated in the synthesis of the actual observational errors. In this paper we follow the 
of test galaxies performed below. An alternative and poten- traditional approach of treating a as a relevant constraint 
tially more attractive approach to deal with the errors is to in the synthesis process, 
marginalize eq. (^) over a. This yields more conservative es- 
timates for the parameters, insofar as they are independent In order to insure a realistic correspondence between 

the values of the errors and the actual quality of the data 
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we have adopted a prescription for the cr's based on our 
experience in dealing with the measurement of W's and 
C's in galaxy spectra (Cid Fernandes, Storchi-Bergmann & 
Schmitt 1998). Two sources of uncertainty affect the mea- 
surements of W's: (1) the noise within the line window, 
and (2) the uncertainty in the positioning of the pseudo- 
continuum fluxes Ck at the pivot wavelengths, which affect 
Wj because they define the continuum over the line. The 
effect of the first of these sources upon o-{Wj) is straight- 
forwardly obtained through standard propagation of errors 
by specifying the noise level at Xj , the spectral resolution SX, 
and the size Aj of the pre-defined window over which Wj 
is measured. In Bica & AUoin (1986a) system the pseudo 
continuum is defined interactively by visual inspection of 
the spectrum around the pivot A's, which makes the es- 
timation of a{Ck) non straight-forward. In Cid Fernandes 
et al. (1998) we found these errors could be roughly scaled 
from the signal-to-noise ratio measured in nearby 'line-free' 
regions, such that the 5'/A'^ in the pseudo continuum is 
(5'/iV)pc ~ ppc{S/N)x, with pre ~ 2-3. This leads to the 
following expression for the a{Wj)'s: 



<j\W,) = A,5x{S/N)l^ + (A, - W,fp^l{S/N)- 



(9) 



where the two terms correspond to the two sources of un- 
certainty discussed above. For simplicity, we further assume 
that the noise in the observed spectrum is constant for all 
A's, which results in {S/N)\ = C\{S/N)5sto- The errors in 
the Cfe's thus become also constant: 



a{Ck)=PphiS/N)^s\o 



(10) 



Besides yielding realistic values for the cr's, this recipe 
has the advantage of quantifying the quality of the data in 
terms of a single quantity: The signal-to-noise ratio at the 
normalization wavelength 5870 A. Of secondary importance 
are the spectral resolution and ppc, which in the simulations 
below are kept fixed at 5 A and 3 respectively. 

In order to provide a quantitative notion of how a 
value of S/N translates into cr(Wj)'s, we have computed 
these quantities for the observables corresponding to spec- 
tral groups SI and S7 of B88, typical of the nuclei of early 
(red) and late (blue) type spiral galaxies respectively. For 
SI our error recipe yields cr(WcaiiK) = 36A/(S'/A'')587o and 
o-(WMgi) ~ IQA./ (S/N) 5S70, so errors of ~ 1 A or less 
are achieved with {S/N)5S70 > 30. For S7 these lines have 
a{W) ~ 15A/(S'/iV) 5870- 



3.3 Exploration of the parameter space: Sampling 
methods 

The exploration of the (x, Ay) space can be approached 
from two alternative perspectives: (1) A minimization prob- 
lem, employing methods to search for the global mini- 
mum (maximum likelihood). Minimization techniques were 
explored in S91 and Pelat (1997) for the synthesis of W's 
only. (2) A sampling problem, where one seeks to map out 
the posterior probability distribution (eq. |^ . 

The methods discussed in this paper are sampling meth- 
ods, which do not explicitly search for a minimum. In the 
limit of small errors, however, one expects the posterior 
maps to peak at the most likely model. For large errors, on 
the other hand, the best model sampled may be well off the 



true global minimum, but in this case the probability distri- 
bution is so broad that the very meaning of 'best model' is 
questionable. In such cases it is arguably more important to 
estimate plausible ranges for the parameters than to carry 
out a refined search for the most likely values. 

Our main motivation to explore the sampling approach 
is that most applications of EPS with Bica's base to date 
focused on the estimation of mean parameters, whose mean- 
ing cannot be directly ascertained without estimating their 
uncertainty and/or biases. A critical re-evaluation of this 
method is thus crucial to establish to which level of detail 
the results of this popular technique can be trusted. 



3.3.1 Uniform sampling 

In order to compute the individual posterior probabilities 
for each parameter, we first note that eq. (p|) for the Xi 
posterior is simply the mean probability P(x, f, 7i) 

over the space spanned by all possible values of Xj^i and 
Av for a fixed Xi. The simplest method to compute such 
probabilities is to divide the (x, Av) space into a uniform 
grid with Aa; steps for the population fractions and AAy 
for the reddening parameter. One then approximates eq. (^ 
by the finite sum 



P{x^\V,a,H) 



-£(x3,Ai,,J 



(11) 



where s denotes a point (a 'state') in the 13-D grid and 
the S{xi — Xi^s) term retains only those grid points where 
population i contributes with fraction Xi of the light. 

This discrete sweeping of the parameter space is easy 
to implement, and was in fact the recipe followed by most 
works with Bica's base to date (see section ^). A serious 
drawback of uniform sampling is its computational price. 
The number of grid cells is an extremely steep function of 
the resolution Ax. A coarse Ax — 10% sampling yields 
3.5 X 10^ x-points, increasing to ~ 8 x 10^ for Ax = 5% 
and more than 4 x 10^^ for a 2% resolution! To obtain the 
grid size one has to further multiply these numbers by the 
number of points in the Ay grid, whose limits have to be 
pre-defined. Schmitt etal. (1999), for instance, performed a 
synthesis study with Aa; — 5% and AAv ~ 0.06 for Ay 
between and 1.5, which yields 2 x 10^ grid points, perhaps 
the finest grid ever used in EPS with this 12 elements base. 
One obviously always aim for the best resolution possible, 
but a priori estimates of what resolution is necessary are not 
straight-forward. A 'good' resolution should produce prob- 
ability distributions that are smooth over scales of Aa; and 
AAy. The width of these distributions is of course a func- 
tion of the observational errors and of the set of observables 
being synthesized. 



3.3.2 Metropolis sampling 

A further drawback of uniform sampling is that the algo- 
rithm spends most of the time in regions of negligible prob- 
ability, which contribute little to the integral in equation (^, 
or its numerical version (eq. ^l|). A more efficient sweep of 
the parameter space would be to use a sampling scheme 
which traces the full probability distribution (eq. One 
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such 'importance sampling' scheme is the Metropolis algo- 
rithm (Metropolis etal. 1953), which samples preferentially 
regions where the posterior is large. This is the parameter- 
space exploration method adopted in this work. 

Our implementation of the Metropolis sampler was as 
follows. Starting from an initial arbitrary point, at each it- 
eration s we pick one of the 13 variables at random and 
change it by an uniform deviate ranging from — e to +e, 
producing a new state s -|- 1. If the variable is one of the 12 
Xi's, the whole set is renormalized to unit sum. Moves to- 
wards unphysical values {xi < or > 1, Av < 0) were trun- 
cated. Downhill moves (i.e., towards smaller x^) are always 
accepted, whilst changes to less likely states are accepted 
with probability e~^^'+'^~^'\ thus avoiding trapping onto 
local minima. This scheme, widely employed in statistical 
mechanics (Press et al. 1992; MacKay 2001 and references 
therein), produces a distribution which tends to the correct 
one as the number of samples Na increase. There is, how- 
ever, no universal prescription to choose optimal values for 
e or Ns, so some experimentation is needed. 

3.4 Performance tests 

In order to test our Metropolis EPS code we performed a 
series of simulations for artificial galaxies generated out of 
the base. Two sets of test galaxies were used. The first is 
composed of the 26 galaxies used by S91 (see their Table 5) 
to test their minimization algorithm. As they did not syn- 
thesize colors, we have used a fixed value of Ay = 2/3 to 
redden all their galaxies. The second set is composed of 100 
test-galaxies whose population vectors and extinctions were 
generated randomly using a scheme which insures that in 
many cases a few components dominate the light. All sim- 
ulations presented in this paper were performed sampling 
Ns = 10* states, with the 'visitation parameter' e set to 
0.005 for both the Xi's and Ay - Experiments were also per- 
formed changing these quantities, and found to give nearly 
identical results. We note that less samples can be used if 
one is only interested in the mean (x, Av), as this naturally 
converges much faster than its full probability distribution. 
Simulations were performed for severals values of {S/N)5S70 
and different sets of observables. Here we restrain the dis- 
cussion to the illustration of the sampling method. 

In Fig. |l| we plot the individual projected posterior 
probability distributions for each of the 12 a^i's plus Av for 
two of S91's test galaxies composed of widely different stel- 
lar populations. The posteriors were computed synthesizing 
all 9 lines and 7 colors in the base ('Set A', as defined be- 
low) and using (S/A^)5g7o ~ 1000 to define the errors in the 
observables. This is clearly an unrealistic, highly idealized 
situation, but it serves as a test-bed for the method in the 
limit of perfect data. The simulations were started from the 
Xi = 1/12 and Av ~ 0.5 point. The plot shows that the 
posterior probabilities are sharp functions peaking at the 
expected input value, marked by upper arrows, demonstrat- 
ing that the Metropolis sampler converges to the most likely 
region. One notes, however, that the strong components in 
the old age bin (such as Xi, X2 and x^ in the 'old galaxy' 
example, drawn as solid lines) have broader individual pos- 
teriors. This effect is further discussed in the next section. 

In the opposite limit of large errors, a cx> (the 'infi- 
nite temperature' regime, as it would be called in statistical 



mechanics) the likelihood loses its discriminating power and 
becomes ^ constant over the whole parameter-space. Eq. (^ 
then tends to P{xi) oc (1 — li)^" for a = 12 base. In this 
regime the data does not constrain the model at all, and 
this expression for P{xi) simply refiects the volume of the 
x-space spanned by a given value of Xi. Not surprisingly, 
this distribution peaks at Xi = 0, and has a mean value 
Ti = 1/12. We have verified that the Metropolis scheme re- 
covers this distribution satisfactorily as S/N 0, which 
serves as a further test of the adequacy of this sampling 
method. The skewing of several of the individual posteriors 
plotted in Figs. ^ and ^ towards small Xi is partly due to 
this effect, which drags the components to an intrinsically 
larger region of x-space. 

The evolution of the Metropolis walk through the 13-D 
(x, Av) space is illustrated in Fig. |^. For clarity, only 6 com- 
ponents of X and just one state for each 10* steps are shown. 
Note that projected P{xi) posteriors such as those in Fig. ^ 
are essentially histograms of the Xi states sampled in walks 
as that in Fig. |^, but weighted by their exp— x(^)^/2 like- 
lihood. The input values of the components shown are indi- 
cated at the right edge of the plot for comparison. This par- 
ticular example corresponds to a simulation with S/N = 60 
and with all 16 observables synthesized. One sees relatively 
large excursions of the parameters from their true values. 
Furthermore, anti-correlations between some of the compo- 
nents (xi and X2, xg, and xg) are clearly visible. This 'mirror 
effect' becomes much more pronounced as S/N increases, 
whereas for low S/N the broadening of the likelihood func- 
tion washes out such correlations. This is a consequence of 
redundancies between several of the base components, as 
further discussed below. The top panels in Fig. |^ show the 
evolution of the a-nd how two of the synthesized observ- 
ables, VK(CallK) and C4020, oscillate very much within their 
respective ±1 sigma error ranges. One therefore anticipates 
an excellent fit to the data, and indeed the obtained for 
the meani-x., Av) 'solution' in this example is just 0.8 (see 
also Fig. 0). 

Overall, we conclude that the Metropolis algorithm is 
an efficient sampler of the parameter space for EPS prob- 
lems. It provides a much more continuous mapping of x and 
Av than would be possible with a uniform grid. Further- 
more, it is substantially more efficient computationally. A 
uniform grid with 10* points would only sample each Xi at 
~ 7% steps and 10 values for Av- By concentrating on the 
relevant regions of the parameter space, we achieve a much 
better resolution, as can be judged by the smoothness of the 
posterior distributions on scales of 1-2%. 

One possible variation of this scheme is to combine the 
uniform grid and Metropolis routines, starting a Markov 
chain from each point on a uniform grid. Experiments with 
this alternative scheme were performed, and found to pro- 
duce essentially the same results. This is to be expected, 
since we verified that the starting point has little effect upon 
the resulting posterior distributions. In fact, we performed 
a series of simulations using the mput parameters as the 
starting point and obtained identical results. Though further 
improvements are possible, the simple Metropolis sampler 
already provides a substantial improvement over previous 
works with this base. 
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Figure 1. Individual posterior probability distributions for each of the components of the population vector x and the extinction, 
computed for an old population dominated test galaxy (solid lines) and a 'younger' galaxy (dotted lines). The top arrows indicate the 
input values of the corresponding parameter in each panel. The numbering of the x components follows the convention in Table 1. All 16 
observables in Bica's base (9 absorption lines + 7 continuum points) were synthesized in this example and a S/N of 1000 was adopted 
to illustrate the convergence of the Metropolis sampler in the limit of (almost) perfect data. 



4 TESTS OF THE BASE 



performed for {S/N)5870 = 10, 30, 60, 
and three different sets of observables: 



100, 300 and 1000, 



Having developed a probabilistic formulation of EPS and a 
new sampling method, we now present a series of numer- 
ical experiments designed to evaluate the actual ability to 
recover the stellar population mix of galaxies using Bica's 
base. The tests were performed with the same series of fic- 
tional galaxies described in section 3.4. Simulations were 



• Set A: All 16 observables in the base. 

• Set B: Only the 9 absorption line equivalent widths. 

• Set C: The W's of Call K, CN, G-band and Mgl-^MgH 
plus the continuum fluxes at 3660, 4020, 4510 and 6630, all 
normalized to 5870 A. 
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Figure 2. Illustration of the 'step by step' Metropolis walk through the parameter space. The six bottom panels illustrate the evolution 
of selected components of the population vector (in %), whose input Xi values are indicated at the right edge of each plot. The values of 
the equivalent width of CallK and of the color C4020 a^re plotted for each step in the third and second panels from the top. The dotted 
horizontal lines in these plots mark the itlcr range around the observed value. The top panel shows the for each sampled state. For 
clarity, only the first few million steps are show, and only 1 for every 10000 sampled states is actually plotted. All 16 observables in the 
base and S/N = 60 were used in this example. 



Set A is the ideal, as it uses all information in the base. 
Set B is used as a test of how much is actually gained by 
synthesizing colors as well as W's, while Set C is composed 
of the observables used by Cid Fernandes etal. (1998) and 
Schmitt et al. (1999) to characterize the stellar content of ac- 
tive galaxies and their hosts through long-slit spectroscopy. 
These data will be used to perform a spatially resolved EPS 
study in a future communication. 



4.1 Effects of the errors in the observables 

We first investigate the effects of the measurement uncer- 
tainties in the observables upon the results of the synthesis. 
Qualitatively, one expects the degradation of the data qual- 
ity to broaden the probability distributions. Furthermore, a 
shift and skewing of the P{xi)'s of intrinsically large compo- 
nents towards small Xi's is also expected for low S/N due to 
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intrinsically larger number of state s in this region of x (the 
'dragging effect' discussed in ^3.4). This is exactly what is 
observed in Fig. ^, where we overplot the individual pos- 
teriors for S/N = 1000, 100 and 10 for one particular test 
galaxy of S91, synthesized with Set A observables. As in 
Fig. nl, the true parameters are indicated by the top arrows. 
The broadening of posteriors in this case can be safely at- 
tributed to the errors alone, since, as demonstrated by Pelat 
(1998) most of the test galaxies in S91 (including the one in 
Fig. y) have a unique algebraic solution in a Vl^-only EPS 
problem, despite the fact that the number of degrees of free- 
dom (11 if we do not model the colors) falls short of the 
number of observables (9 W's). 

Despite the large number of constraints in Set A, one 
sees that the individual posterior distributions are substan- 
tially broad even for a S/N = 100 spectrum. The average 
standard deviation of the Xi's for this combination of S/N 
and observables is 3% over all simulations, but can reach 
more than 10% in individual components. Furthermore, the 
dominant components, 2, 6 and 9 in this example, are all 
badly affected by the errors. Fig. ^ shows that decreasing 
data quality progressively shifts P{x2) towards smaller val- 
ues, whereas the opposite happens with xi . This compensa- 
tion is also seen between components in the 1 Gyr (15, X(, 
and 2:7), 100 Myr (xg and 19), and 10 Myr (xio and a;ii) 
age bins and is the same 'mirror' effect observed before in 
Fig. ^. A compelling visualization of these compensations 
is given in Fig. ^, where we plot the individual Metropolis 
states for the test galaxy employed in Fig. ^ in a a;^ x Xj 
matrix. The S/N of 300 in this plot was chosen for clar- 
ity purposes; the same structure is present for other values, 
with the scatter around the strongest correlations increasing 
steadily as S/N decreases, to the point that at S/N — 10 
the anti-correlations between components 1 and 2, or 8 and 
9 become clouds similar to xi x xi, in Fig. ^. All anticorre- 
lations identified in Fig. |^ {xi x X2, xe x xs, etc.) are clearly 
represented in this graphical version of the correlation ma- 
trix. 

These effects are rooted in the internal structure of the 
base. By definition, the base elements must be linearly in- 
dependent, and the S91 base complies with this algebraic 
condition. Yet, some of its elements are practically linearly 
dependent, in the sense that combinations of other elements 
can recover them to a high degree of accuracy. As the noise 
increases, the residuals between representing an element by 
its exact M^'s and C's and a combination including other ele- 
ments with more extreme values for the observables become 
statistically insignificant, explaining the compensations de- 
tected above. To demonstrate this we synthesized the base 
elements themselves with our code. Results for Set A ob- 
servables and three different values of S/N are shown in 
Fig. 1^. Each panel corresponds to one of these 12 extreme 
test galaxies (labeled B1...B12), with the mean synthetic 
values of each of the components plotted along the vertical 
axis. The empty, shaded and filled histograms correspond to 
S/N = 100, 30 and 10 respectively. Some of the components, 
most notably numbers 2, 3 and particularly 6, are synthe- 
sized with large contributions of others even for S/N = 100. 
For the xe = 100% model (top left panel of Fig. q), for in- 
stance, we find xE = 60% for this high S/N, with 35 of the 
remaining 40% redistributed among components 5 and 7. 
Yet, the 16 observables are very well reproduced by this 



combination, with a total just 3.8. This spreading 

of the light fractions becomes more pronounced for lower 
S/N's, to the point that at S/N = 10, components like 8 
and 9, or 10 and 11 cannot be distinguished anymore, and 
end up dividing their contributions in roughly equal shares, 
with residuals spread over other components. The spread is 
much smaller in the synthesis of components 1, 4, 5 and 7, 
which represent extreme metallicities within age groups. For 
intermediate Z systems (models B2, B3 and B6), however, 
these extreme components attract much of the percentage 
light contribution spilled over from neighboring populations 
of same age, so that their individual contributions in a real 
stellar population mix are as unreliable as the others. 

This experiment demonstrates that for practical pur- 
poses the base is still linearly dependent, at least in a statis- 
tical sense, despite the effort of S91 in reducing the highly 
redundant original base of B88. This 'noise- induced', or 'sta- 
tistical linear dependence', as we may call it, could indeed 
be inferred from the PCA analysis of S91, which showed 
that very little of the variance is associated with the last 
few eigenvectors. Within our probabilistic formulation, this 
'statistical linear dependence' defines the structure of the 
likelihood function, which spreads first along directions in 
x-space producing similar observables, carrying the mean x 
to more densely populated regions. This effect explains the 
redistribution of the probability observed in Figs. ^, ^ and 
all other simulations. 

It is therefore extremely difficult to retrieve accurately 
all the input stellar population parameters even for excellent 
data. This is further illustrated in Fig. ^ where the input 
parameters are plotted against their mean values, as derived 
from the Metropolis runs. One hundred artificial galaxies are 
plotted in each panel, with open circles corresponding to Set 
A observables and S/N = 1000, and crosses indicating the 
results for Set A and S/N = 60. A systematic underestima- 
tion of strong components is seen for 5'/A*' = 60, a good spec- 
trum by any standard. This happens at the expense of an 
overestimation of the weaker components, which produces 
the large scatter seen in the bottom left corners of Fig. ^ 
The average uncertainties in the Xi's are of order a{xi) ~ 3- 
4% for S/N = 60. These are not enough to account for the 
large differences between input and output Xi^s depicted in 
Fig. so the errors induce a true bias in x. Note also that, 
in accordance with the discussion above, this bias is smaller 
for models with strong contributions of components 1, 4, 5 
and 12, because of their extreme locations in the space of 
observables. 

Only for unrealistically large S/N's, which far exceed 
the quality of the data used to build the base in the first 
place, one is able to break the 'statistical linear dependence' 
of the base by distinguishing fine details in the observables. 
This should not be interpreted as a failure of the method, 
as what the code actually synthesizes are the observables! 
These are very precisely reproduced by the mean (x. Ay), as 
illustrated in Fig. gfor the W's of Call K, CN, G-band and 
Mgl and two continuum colors. The remaining observables, 
not show in the figure, are equally well reproduced. Further- 
more, whilst S/N > 300 is needed to recover accurately the 
model parameters from the mean solution, the agreement 
between the synthetic and measured observables is excellent 
for any S/N. Obviously, an even better agreement is ob- 
tained if instead of {x.,Av) the best model sampled during 
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Set A & S/N = 1000, 100, 10 
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Figure 3. Effect of SjN upon the individual probability distributions of the synthesis parameters for Set A observables. The different lines 
correspond to SjN = 1000 (solid), 100 (dotted) and 10 (dashed). Notice how the probability of x components with a large contribution 
(like X2, x% and in this example) are progressively redistributed among other components as the noise increases. 



the Metropolis excursion is used to reconstruct the observ- 
ables, as illustrated in the fourth column of plots in Fig. ^ 
(see §4.5). 



In conclusion, this comparison of input and output ob- 
servables shows that the spreading of the probability and 
the confusion between components seen in Figs. I and | IS 
not an artifact of the method, but a consequence of the in- 
ternal structure of the spectral base. As found in other EPS 
studies (O'Connell 1996 and references therein) synthesiz- 
ing the observables is one thing, but trusting the detailed 



star-formation history and chemical evolution implied by 
the synthesis is an altogether different story. On the pos- 
itive side, a notable fact about Fig. ^ is that, consulting 
Table |^, one realizes that the reshuffling of the strength of 
the components occurs preferentially among populations of 
same age, an effect also clearly seen in Fig. El Though some 
redistribution among base elements of different age and Z 
in the older 1 and 10 Gyr bins due to the age-Z degeneracy 
also occurs, this indicates that t he a ge distribution may be 
well recovered by the synthesis (§4.3). 
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Figure 4. 'Correlation Matrix' for the test galaxy in Fig. g| Each Xi X Xj panel shows one every 10000 of the first lO''^ steps of a Metropolis 
run for {S/N)ss7o = 300 and Set A observables. The values at the bottom of each row indicate the input value of the corresponding 
component (xi = 0, X2 = 20%, . . .xn = 10%, xi2 = 10%). These are also indicated by the solid circles in each panel. Small and big tick 
marks are spaced by 2 and 10% respectively in all plots. Note the strong anti-correlations between adjacent components of same age. 



4.2 Effects of different sets of observables 

The effects of which set of observables is used in the synthesis 
are illustrated in Fig. |^, where results for Sets A, B and C 
for the same test galaxy as in Fig. ^ and S/N = 1000 are 
compared. The figure shows that it is very advantageous to 
synthesize colors along with equivalent widths, despite the 
fact that one increases the dimensionality of the problem 
with the inclusion of Ay - This can be seen by comparing the 
performances of Sets A and B. The information contained 



in the colors not only improves the estimation of x but also 
allows a very good determination of Av, which, unlike the 
population vector, is always well recovered by the synthesis. 

Fig. ^also shows how important it is to provide as many 
observational constraints as possible to constrain the syn- 
thesis, as Set A recovers the parameters much better than 
Sets B and C. Decreasing the number of observables in the 
synthesis thus have the same overall effect as decreasing the 
S/N. The width of the posteriors for Set C are partially 
attributed to the algebraic degeneracy in this set (8 observ- 
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Figure 5. Synthesis of tlie base elements. Each panel shows the mean synthetic population vector (vertical axis) for 12 test galaxies 
whose 16 observables were generated with Xi = 100%, i = 1, ... 12 (labeled Bl. . . B12 respectively). Empty, shaded and filled histograms 
correspond to S/N = 100, 30 and 10 respectively. Ideally, all histograms should be concentrated in the component used to generate the 
observables. In practice, some bases elements, most notably 2, 3 and 6, are well synthesized by combinations of others even for large 
S/N. 



ables and 12 degrees of freedom), which implies the exis- 
tence of exact solutions in a sub-space of (x,^^). That is 
not the case of Set A for the test galaxy in this example 
(Pelat 1998), but there, as for other sets, the 'statistical lin- 
ear dependence' of the base is very efficient in spreading the 
likelihood even for such an idealized S/N ratio. The excellent 
agreement between the 'observed' and synthesized observ- 
ables even for much worse S/N (Fig. ^) illustrates this point. 
Algebraic degeneracy is not as critical as the 'statistical de- 



generacy' induced by the combined effects of noise, limited 
data and the internal correlations in the base. 

Studies which synthesize only W^'s sometimes use the 
resulting population vector to compute a predicted contin- 
uum shape and derive Ay a posteriori through the compar- 
ison with the observed continuum (e.g., B88; Boisson etal. 
2000). From our Set B simulations, we find that this pro- 
cedure yields values of Av 2-3 times less accurate than ob- 
tained with Set A. For S/N — 60, for instance. Set A yields 
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Figure 6. Input X output plots for all 13 individual parameters in the synthesis. Open circles map the results for Set A and 3/N = 1000, 
while crosses correspond to Set A and S/N = 60. The output parameters are the mean values sampled from the global likelihood function. 
The dotted, diagonal lines in each panel mark the 'input = output' line. Error bars on the model parameters are not shown for clarity. 
Notice the systematic underestimation of large x^'s and the corresponding overestimation of weak proportions for S/N = 60. 



an rms dispersion in Ay of 0.06 mag around the true value, 
whereas the dispersion for set B is 0.13 mag. Since this is 
a relatively small difference, these experiments validate the 
a posteriori computation of the extinction. Colors are more 
useful to constrain the population vector than to estimate 

Ay. 



4.3 Age grouped results 

The results of the tests above reveal a striking difficulty to 
accurately determine all 12 base components of Bica's base 
in the presence of even very modest noise and/or when not 
all base observables are available for the synthesis. As in- 
dicated by the simulations, both the measurement errors 
and the use of reduced sets of observables act primarily in 
the sense of spreading a strong contribution in one compo- 
nent preferentially among base elements of same age. Group- 
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Synthetic Quantities 



Figure 7. Input X output observables for a series of 100 test-galaxies synthesized with the Metropolis-EPS sampler. The three left 
columns of plots correspond to the three indicated combinations of the set of observables and S/N. For these plots the x-axis quantities 
were computed with the mean (x,Av) solution. The plots on the right column were made with Set C and S/N = 30 (as those in the 
third column), but the synthetic observables in this cases were computed with the best model found during the Metropolis runs. Only 
the equivalent widths of four absorption features and two colors exhibited, but results for the remaining observables are equally good. 
The vertical and horizontal scales are the same for each panel, and the input = output is indicated by a dotted line. Notice the excellent 
agreement, despite the fact that the input parameters are not accurately recovered by the synthesis process. 



ing the population vector in age bins should thus produce 
more robust results. This expectation is confirmed in Fig. ^ 
where, analogous to what was done in Fig. |^, we plot the 
input values of the Xi's against the output afl's, but now for 
the five age binned groups, also for Set A simulations. The 
left panels correspond to the same values of S/N used in 
Fig. H, so that one can appreciate the enormous improve- 



ment achieved by grouping equal age elements. The right 
panels show that good results are also obtained for S/N of 
30 (circles), and that young populations (< 100 Myr) are 
reasonably well traced even for S/N as low as 10 (crosses). 
Though some deviations survive, specially for the smaller 
S/N's and the older age bins, the age-binned synthesis re- 
sults are much more reliable than the results for the indi- 
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Figure 8. Analogous to Fig. but exploring the effect of different sets of observables upon the synthesis parameters. Solid lines 
correspond to Set A (all observables), dotted lines to Set B (only 9 equivalent widths) and dashed lines to Set C (4 lines and 4 colors). 



vidual components, and, most importantly, do not require 
unrealistic 5'/A'^. 

This conclusion also holds for Sets B and C, as illus- 
trated in Fig. Naturally, the restricted input data in these 
sets translates into a loss of information, and the agreement 
is not as good as for Set A. Note that Set C does a much 
better job in recovering populations of 100 Myr or less than 
Set B, a further example of how useful colors are in the syn- 
thesis. This happens because the color information (present 
in Set C but absent in B) impose stronger constraints on 
young populations than in older systems. 



The reliability of age-grouped proportions contrasts 
with the badly constrained and biased results for the in- 
dividual components of x, and is one of the main results of 
our study. This result puts in perspective all previous inter- 
pretations of the synthesis with this base. This fundamental 
limitation of the base was in fact previously known, but was 
never studied in detail. Schmitt etal. (1999), for instance, 
preferred to summarize the results of their synthesis study 
in terms of age grouped proportions for populations younger 
than 1 Gyr, in consonance with the results above. Still, in 
that work we kept a distinction between the different Z's in 
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Figure 9. Age-binned light fractions in 100 input test ealaxies versus the corresponding synthetic means. Left panels: Circles correspond 
to S/N = 60, and crosses to S/N = 1000, as in Fig. fel Right panels: Circles correspond to S/N = 30, and crosses to S/N = 10. All 
results are for Set A observables. 



the 10 Gyr age-group. At this level of detail the population 
fractions are very uncertain. In fact, since Set C observables 
and S/N ~ 30-60 were used in that study. Fig. ^ indicates 
that it would be safer to group their 1 and 10 Gyr pro- 
portions. Yet, the proportions for the younger (< 100 Myr) 
populations found by Schmitt etal. (1999), and which con- 
stituted the focus of that paper, are trustworthy (Fig. Fui). 



4-3.1 The age x metalUcity degeneracy 

A close inspection of the S/N = 10 and 30 results in Fig. ^ 
reveals a compensation effect between the 1 and 10 Gyr age 
bins, albeit at a much smaller level than for the individual 
components (Fig. This effect is much more pronounced 
in Sets B and C, as seen in Fig. |l^. 

The confusion between the 1 and 10 Gyr populations 
is related to the a.ge-Z degeneracy, which sets in precisely 
in this age range. This well known effect (O'Connell 1986, 
1994; Worthey 1994) is also present in Bica's base and affects 
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Figure 10. Input X output age binned x fractions for Set B (left) and Set C (right) observables. Circles correspond to S/N = 100, and 
crosses to S/N = 30. 



the synthesis results. To visualize and quantify this effect 
we have computed the average log Z and log age for both 
the input and synthetic x for our test galaxies. Given that 
X is a light fraction at 5870 A, these averages ultimately 
represent flux-weighted ages and metallicities, which serve 
our current illustrative purposes, despite their ambiguous 
physical significance. 

The 'output - input' log age and log^ residuals are 
plotted against each other in Fig. In these diagrams 
A log age > (overestimated age) implies redder colors and 
stronger metal lines, which tend to be counter-balanced by 



the bluer colors and weaker lines resulting from A log Z < 
(underestimated Z). The a,ge-Z degeneracy is most obvious 
in the results for Sets B and C, for which the biases in age 
and Z persist even for S/N = 100. For Set A, on the other 
hand, the residuals are small: ^ 0.2 dex for S/N ^ 60, 
which shows that the age-Z degeneracy can be broken, at 
coarse the level of age and Z resolution offered by the base, 
with enough information and good spectra. Note that this 
remark applies to the global averages defined above, not the 
detailed component by component ages and Z's, which we 
concluded to be highly uncertain (Fig. 
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Figure 11. The age-metallicity degeneracy. The x-axis maps the difference between the output and input log age, with the corresponding 
difference in logZ plotted along the y-axis. Circles correspond to age and Z computed from the mean {x,Av), whereas dots correspond 
to the best model sampled. All 126 test galaxies are plotted. The different panels correspond to combinations of the data quality {S/N, 
labeled in the bottom-left corners) and the set of observables (top-right) used in the synthesis. 



4.4 Metallicity grouped proportions 

The results for Z binned proportions, plotted in Fig. |l2| 
are not nearly as good as those for the age binned groups. 
With the exception of the Z = Zq bin, which is actu- 
ally represented by just one element (Table 1), the scatter 
in the input x output Z-binned proportions is so large for 
all other Z's that one is forced to conclude that no accu- 
rate description of the chemical history of galaxies can be 
afforded with this spectral base. Given this limited Z 'res- 



olution power', one might consider removing intermediate 
components such as 2, 3 and 6. However, as we concluded 
that only the age distribution can be assessed with the base, 
there is no obvious advantage in this further reduction. This 
situation can be improved by providing a priori constraints 
on the occupancy of the age-Z plane, like those imposed by 
B88 based on chemical evolution arguments, since these ef- 
fectively reduce the dimensionality of the base to typically 8 
components, thus alleviating the confusion between compo- 
nents and producing better focused results. The validity of 
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Figure 12. Metallicity binned light fractions in 100 input test 
galcixies versus the corresponding synthetic means. Circles corre- 
spond to S/N = 100, and crosses to S/N = 30, both for Set A 
observables. The 'resolution power' of the base is much better for 
ages (Fig. ^) than for metallicities. 



this prior input external to the synthesis process, nonethe- 
less, has to be evaluated in an object by object basis. In 
this study we followed S91 in not imposing any such extra 
information, as this approach encompasses a larger class of 
evolutionary scenarios and thus contemplates more possible 
applications. 



4.5 Best X Mean Parameters 



mass' of acceptable solutions. The severe biases identified 
in the mean synthetic population vector, however, raise the 
question of whether alternative approaches should be pur- 
sued instead. 

By definition, minimization oriented techniques would 
do a better job at recovering the input parameters. Indeed, 
we verified that the best model population vector sampled in 
the Markov chain is much closer to Xtrue than its expected 
value, an agreement which can only be improved given that 
the Metropolis sampler is not an optimal minimization pro- 
cedure. Though the simulations confirmed our expectation 
that for small errors the mean and best solutions should both 
converge to the true parameters, it was surprising to realize 
how easily x moves away from xtme, and consequently also 
from xtest. Very small uncertainties in the observables, at 
the levels of S/N > 300, are enough to set powerful com- 
pensation effects into operation (e.g.. Figs. ^ and^, shifting 
the components of x away from their true values due to the 
asymmetric redistribution of the probability in x-space. 

Minimization procedures can in principle overcome this 
bias, as illustrated by the success of the S91 and Pelat (1997) 
tests as well as by our (not optimized) results for x^est. A 
further bias in mean solutions is that they always produce 
aJT > 0. Real galaxies, even if synthesizable by the base, are 
likely to have observables outside synthetic domain due to 
measurement errors. As shown by Pelat (1997, 1998), the 
best model in this case lies at the surface of the synthetic 
domain, where very often several of the x-components are 
0. In these aspects, minimization algorithms are more suit- 
able than the more traditional sampling approach. However, 
the high susceptibility of the synthesis to the data quality 
indicates that minor perturbations in the input observables, 
within their errors, may seriously affect the estimation of 
the best model parameters. The stability of the best solu- 
tion is thus likely to be more fragile than that of the mean 
solution, which takes into account the effects of the errors 
in the observables in the trade-off between an optimal and 
a statistically acceptable fit of the data. We have in fact 
detected this effect in a series of simulations, but a more 
careful study, with specialized minimization techniques is 
required to quantify it properly. It would be particularly in- 
teresting to carry out this test with the method developed 
by Pelat (1997) applied to the base used in this work. His 
Monte Carlo tests with a n^. = 10 stellar base and nw = 19 
equivalent widths produced essentially unbiased population 
fractions even for S/N as low as 10. Extending his formalism 
to synthesize galaxy colors would also be desirable. 

In any case, because of the nature of the EPS problem 
and the structure of Bica's base, relatively large discrepan- 
cies in the parameters do not necessarily translate into large 
differences in the synthesized quantities. In this respect, the 
difference between mean and best models is largely aca- 
demic, as the mean solution already does an excellent job in 
fitting the observables (Fig. m. 



Throughout these experiments, we have consistently used 
the mean parameters as an estimate of the result of the 
synthesis. This convention was deliberately chosen to follow 
the more widely employed version of EPS with this base 
and thus to allow a reassessment of previous results. Fur- 
thermore, the mean is a convenient and reasonable way to 
summarize the results insofar as it represents a 'center of 



5 SUMMARY 

In this paper we have revisited the method of Empirical 
Population Synthesis with the aims of improving it both 
at the formal and computational levels, and exploring its 
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efficacy as a tool to probe the population mixture of gala^xies. 
Our results can be divided into three parts. 

(1) A simple probabilistic formulation of the problem 
was presented, which puts this method, traditionally em- 
ployed in a loss foruial way. onto a mathematical footing. It 
was shown how former applications of the method fit into 
the probabilistic formulation, thus providing some a poste- 
riori justification for previous results. 

(2) An importance sampling scheme, based on the 
Metropolis algorithm, was developed and tested. It provides 
a more efficient and smooth mapping of the probability dis- 
tribution of the parameters than is possible with the com- 
monly used uniform grid sampling. 

(3) In the third part of the paper, we applied the for- 
mulation and sampling method to a series of test galax- 
ies constructed out the n* = 12 star-clusters base of S91. 
The tests explored the ability to reconstruct the model pa- 
rameters, which carry a record of the star formation and 
chemical evolution in galEixies, in the presence of (i) obser- 
vational errors and (ii) limited data. This study centered on 
the comparison of input parameters with the mean synthetic 
solution (x, j4v), as this is the most common form of EPS 
in the literature with this base, and a systematic evaluation 
of the consistency of this method has not been carried out 
before. The main results of this study can be summarized 
as follows: 

(a) The allowance for errors in the observables sparkle 
linear dependences within the base elements. This induces 
systematic biases in x, redistributing the probability in com- 
ponents with a large fractional contribution to the integrated 
spectrum among similar components. This happens prefer- 
entially for components of same age. Though this was a pre- 
dictable result, it was surprising to realize how powerful this 
effect is. For any realistic S/N, the bias is such that in most 
cases one cannot trust the estimates of all 12 mdividual com- 
ponents of X to any useful level of confidence. 

(b) Reducing the number of observational constraints, 
due to incomplete coverage of the 3300-9000 A spectral in- 
terval spanned by the base, has the same overall effect as 
reducing the data quality, that is, to broaden, shift and skew 
the probability distribution of the Xi's. We find that synthe- 
sizing colors as well as equivalent widths yields substantially 
better results than synthesizing W's only. The need to ac- 
count for the extinction as an extra parameter is largely 
compensated by the better focused x obtained by consider- 
ing the color information, particularly for populations of 100 
Myr and younger. Unlike for the individual population frac- 
tions, we do not detect systematic biases in the estimates of 
Av. 

(c) Despite the difficulty to retrieve accurately the de- 
tailed population vector, the observables are very precisely 
reproduced by the mean solution. The mapping between the 
parameters and observables spaces is thus highly 'degener- 
ate'. The nature of this 'non-uniqueness' is not algebraic, as 
it happens also when the number of free parameters exceeds 
the numbers of observables and in undetermined cases where 
the solution is known to be unique. Instead, it is related 
to the possibility of mimicking to a statistically indistin- 
guishable level the effects of certain base elements by linear 
combinations of other elements. This happens even for tiny 



observational errors, inducing a 'spill over' of the likelihood 
among base components, carrying the mean solution along. 

(d) The structure of the base reflects the fact that evo- 
lution is the main source of variance in stellar populations. 
Accordingly, we found that grouping the x components in 
their five age bins produces much more reliable results, and 
we recommend this procedure when evaluating the detailed 
population vector obtained in previously published results 
with this base. Some confusion between the 1 and 10 Gyr 
populations also takes place due to the &gc-Z degeneracy, 
but this can be minimized with an adequate spectral cover- 
age (measuring all observables in the base) and S/N ^ 30. 
Grouping the results for the four different Z's in the base 
is also recommended, but caution must be exerted when in- 
terpreting the results since these are much less precise than 
the age-binned proportions. 

The difficulties to retrieve accurately the true stellar 
population parameters stem from two sources: (i) The con- 
vention to represent the results of the synthesis by the mean 
solution, and (ii) limitations imposed noise-induced linear 
dependences in the base. Minimization techniques may be 
instrumental in overcoming the biases in the mean popula- 
tion vector, which ultimately limit the resolution of the 12-D 
base to its 5 age bins. Within our probabilistic formulation, 
other priors can be investigated. Maximum Entropy, for ex- 
ample, is a prior that has been applied successfully to inverse 
problems formally similar to population synthesis, like im- 
age restoration. However, more than improvements on the 
method, which are possible and desirable, further work on 
the base itself will certainly be needed. Indeed, the limited 
age and Z 'resolution power' of the base employed here lies 
at the heart of all difficulties identified in this paper. In- 
cluding other spectral indices, modeling the full spectrum 
and incorporating other relevant information such as mass 
to light ratios in the synthesis process are possibles ways to 
progress in this direction. 
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